Motivation

As far I know, there is no package for image registration in R and no wrapper to open.CV, so I’ll try to write some simple functions to detect keypoints in images, calculate descriptors , find matches and adjust a homomorphic transformation between images.

There are some packages in R for image manipulation and after some test I select “imager” , based on the CImg C++, fast and providing several image processing tools.

Reading and plotting images

Let’s read an image using the “jpeg” package, and transform to a CImg object for manipulation with “imager”

set=10;id=4
nim=paste0(path,"set",set,"_",id,".jpeg")
im=readJPEG(nim)
str(im)
##  num [1:2329, 1:3100, 1:3] 0.729 0.725 0.725 0.718 0.714 ...

the im object is an 2-d numerical array with depth 3 (for the RGB channels) We can work with the full color image or with a channel. The plot command is able to display a CImg object

cim=imresize(as.cimg(im),scale = scl)
## Warning in as.cimg.array(im): Assuming third dimension corresponds to
## colour
str(cim)
##  cimg [1:698, 1:930, 1, 1:3] 0.722 0.712 0.702 0.679 0.71 ...
cimr=imresize(as.cimg(im[,,1]),scale=scl)#Red Channel
str(cimr)
##  cimg [1:698, 1:930, 1, 1] 0.722 0.712 0.702 0.679 0.71 ...
par(mfrow=c(1,2))
plot(cim)
plot(cimr)

Keypoint detection

The first step in the registration process consist in the keypoint detection. The points selected must be robust to translations, rotations and scaling and stables to noise. There is a handful of methodologies, some of them subject to patents like SIFT, but simpler ideas can also work. In this case we can try with a “corner detector”. The Harris detector is based in gradients and gaussian filtering, and is easily implemented with imager. The value of the Harris transformation is high only in well defined corners on the image. We select as keypoints the center of connected regions with Harris over a certain threshod.

#HARRIS - Harris corner detector
Harris<-function(im,sigma=2){
  eps=1.e-10
  ix=imgradient(im,"x")
  iy=imgradient(im,"y")
  ix2=isoblur(ix*ix,sigma,gaussian = T)
  iy2=isoblur(iy*iy,sigma,gaussian = T)
  ixy=isoblur(ix*iy,sigma,gaussian = T)
  (ix2*iy2-ixy*ixy)/(ix2+iy2+eps)
}
cim_Harris=Harris(cimr,sigma=3*scl)


#Detect Keypoints
get.centers <- function(im,thr="99%",sigma=3*scl,bord=30*scl){
  dt <- Harris(im,sigma) %>% imager::threshold(thr) %>% label
  as.data.frame(dt) %>% subset(value>0 ) %>% dplyr::group_by(value) %>% dplyr::summarise(mx=round(mean(x)),my=round(mean(y))) %>% subset(mx>bord & mx<width(im)-bord & my>bord & my<height(im)-bord)
}
par(mfrow=c(1,2))
plot(log(cim_Harris+.0001))
plot(cimr)
# Detect keypoints 
kp=as.data.frame(cimr %>% get.centers(sigma=3*scl,"98%"))[,2:3]
kp %$% points(mx,my,col="red")

head(kp)
##    mx my
## 1 474 69
## 2 509 17
## 3 404 16
## 4 262 10
## 5 328 10
## 6 476 11

Image Orientation

The next step is to build point descriptors. A simple solution is to use a normalized patch around a blured region of each key point. For instance a grid of 9x9 pixels taken from a 20x20 patch, every 5 pixels

stencil <- expand.grid(dx=seq(-20,20,5)*scl,dy=seq(-20,20,5)*scl)

Usually the patch must be oriented using the main patch direction, but in this case it seems to work better to use a global orientation for all patches. To do that we must calculate a global direction for the whole image. This is done using the histogram of gradient (HoG) orientations and selecting the maximum (correcting for +/- 90º rotations)

sigma_b=6*scl
ima_bl=isoblur(cimr,sigma_b,gaussian = T)

pi=3.141592653

ix=imgradient(ima_bl,"x")
iy=imgradient(ima_bl,"y")
ita=atan(iy/ix)*180/pi
hist(ita,breaks=180)

The peak in the histogram of gradients define the main image direction. Some images can have more than one local maxima in the Hog, defining several possible image orientation Let’s define a function to calculate it.

#Return 1 or 2 main global image orientations restricted to (-45,45)

get_orientations<-function(im){
  ix=imgradient(im,"x")
  iy=imgradient(im,"y")
  ita=atan(iy/ix)*180/pi
  iga=table(sample(round(ita*2)/2,200000))
  #plot(iga)
  ma1=max(iga)[1]
  m1=which(iga==ma1)
  theta_1=(as.numeric(names(m1)))
  iga[max((m1-20),0):min((m1+20),length(iga))]=0
  #plot(iga)
  ma2=max(iga)[1]
  m2=which(iga==ma2)
  theta_2=(as.numeric(names(m2)))
  if(theta_1>45) theta_1=theta_1-90
  if(theta_1<(-45))theta_1=theta_1+90
  if(theta_2>45) theta_2=theta_2-90
  if(theta_2<(-45))theta_2=theta_2+90
  if(abs(theta_1-theta_2)>5){
    return(c(theta_1,theta_2))
  }
  else{
    return(theta_1)
  }
}

th=get_orientations(ima_bl)
print(paste0("Global main directions ",th,"º"))
## [1] "Global main directions -9º" "Global main directions 1º"
par(mfrow=c(1,length(th)))
for(theta in th)plot(imrotate(cimr,-theta))

Keypoint descriptor

Now we can extract the descriptor from a rotated patch, asuring all path are equally oriented across images. If the image had more than one orientation, we have to calculate one patch for every possible orientation. Lets plot some patches

stencil <- expand.grid(dx=round(seq(-20,20,5)*scl),dy=round(seq(-20,20,5)*scl))
stencil_ext <- expand.grid(dx=round(seq(-30*scl,30*scl,1)),dy=round(seq(-30*scl,30*scl,1)))
par(mfrow=c(3,3))


for(i in c(12,50,100)){ #3 random keypoints
  pm=get.stencil(cimr,stencil_ext,x=kp[i,1],y=kp[i,2])
  pm=as.cimg(pm)
  plot(pm)
  imr=imrotate(pm,-th[1])
  plot(imr)
  ww=round(width(imr)/2)
  desc=get.stencil(imr,stencil,x=ww,y=ww)
  plot(as.cimg(desc))
}

#Get oriented descriptors
get_descriptor_oriented<-function(im,theta,v){
  pm=get.stencil(im,stencil_ext,x=v[,1],y=v[,2])
  w=sqrt(length(pm))
  pm=as.cimg(pm,x=w,y=w)
  imr=imrotate(pm,-theta)
  ww=round(width(imr)/2)
  get.stencil(imr,stencil,x=ww,y=ww)
}

Matching of keypoints

Let’s try now with a pair of images:

ida=1;idb=4
nim=paste0(path,"set",set,"_",ida,".jpeg")
ima=imresize(as.cimg(readJPEG(nim)[,1:3099,1]),scale=scl)
nim=paste0(path,"set",set,"_",idb,".jpeg")
imb=imresize(as.cimg(readJPEG(nim)[,1:3099,1]),scale=scl)

sigma_b=6*scl

ima_bl=isoblur(ima,sigma_b,gaussian = T)
imb_bl=isoblur(imb,sigma_b,gaussian = T)


tha=get_orientations(ima_bl)
thb=get_orientations(imb_bl)

par(mfrow=c(1,2))

plot(imrotate(ima,-tha[1]))
plot(imrotate(imb,-thb[1]))

par(mfrow=c(1,2))

plot(ima)
# Detect keypoints 
kpa=as.data.frame(ima %>% get.centers(sigma=3*scl,"98%"))[,2:3]
kpa %$% points(mx,my,col="red")

plot(imb)
kpb=as.data.frame(imb %>% get.centers(sigma=3*scl,"98%"))[,2:3]
kpb %$% points(mx,my,col="red")

feata=NULL
for(theta in tha){
  dfa<-alply(kpa,1,function(v){ ss=get_descriptor_oriented(ima_bl,theta,v)}) %>% do.call(rbind,.)
  dfa=as.data.frame(t(apply(dfa,1,scale)))
  feata <- rbind(feata,dfa)
}
featb=NULL
for(theta in thb){
  dfb<- alply(kpb,1,function(v){ ss=get_descriptor_oriented(imb_bl,theta,v)})  %>% do.call(rbind,.)
  dfb=as.data.frame(t(apply(dfb,1,scale)))
  featb <- rbind(featb,dfb)
}

Now feata and featb contains the descriptors for both images. The next step is to find approximate matching between them. We will use knn to find points close in the feature space: we’ll kept first neighbors when the second neighbor distance is bigger within a factor

require(FNN)

kk<-get.knnx(data=feata, query=featb, k=2, algorithm ="kd_tree" )
if(length(thb)==1){
  lpb=c(1:nrow(kpb))
}else{
  lpb=c(c(1:nrow(kpb),c(1:nrow(kpb))))
}
if(length(tha)==2)kpa=rbind(kpa,kpa)

mask=(kk$nn.dist[,1]/kk$nn.dist[,2]<.8)
match=cbind(kk$nn.index[mask,1],lpb[mask])               

p1=as.matrix(kpa[match[,1],])
p2=as.matrix(kpb[match[,2],])

par(mfrow=c(1,1))
plot(kk$nn.dist[,1],kk$nn.dis[,2],pch='.')
points(kk$nn.dist[mask,1],kk$nn.dis[mask,2],pch='o',col="red")

RANSAC (Random Sample Consensus)

Probably, most of the matchings are spurious, so we need a robust method to fit the affine transformation from one set of points to the other. With RANSAC, we select ramdomly sets of 4 points for adjusting an homography, and record the number of the remaining points that agree with this parameters within an error. Then we select as “inliers” the most populated set and fit a final homography only to them.

Some function definitions:

# Estimate a homography h from points in P to points in p
est_homograph<-function(P,p){
  n=nrow(P)
  hh=NULL
  for(i in 1:n){
    a=t(c(p[i,],1))
    b=t(c(0,0,0))
    c=P[i,]
    d=-c%*%a
    hh=rbind(hh,cbind(rbind(c(a,b),c(b,a)),d))
  }
  h=t(matrix(svd(hh,nv=ncol(hh))$v[,9],nrow=3,ncol=3))
}

#Apply homographyh to points in p
apply_homograph<-function(h,p){
  p1=t(cbind(p,1))
  q1=t(h%*%p1)
  q1=q1/q1[,3]
  q1[,1:2]
}

#Robust homography estimation from p1 to p2. Return h and the list of inliers
ransac<-function(p1,p2,thresh=100,N=1000){
  n=nrow(p1)
  set.seed(12345)
  sn=c(1:n)
  flag=matrix(0,nrow=N,ncol=n)
  for(i in 1:N){
    smpl=sample(sn,4)
    pp1=p1[smpl,]
    pp2=p2[smpl,]
    h=est_homograph(pp2,pp1)
    p=apply_homograph(h,p1)
    d=rowSums((p-p2)^2)
    flag[i,]=as.numeric(d<thresh)
  }
  sinliers=rowSums(flag)
  sinliers=sinliers[!is.na(sinliers)]
  imax=which(sinliers==max(sinliers))[1]
  inliers=sn[flag[imax,]==1]
  h=est_homograph(p2[inliers,],p1[inliers,])
  list(h,inliers)
}

Now, we calculate the homography, and show the inliers in green

hh=ransac(p1[,1:2],p2[,1:2],100,5000)

h=hh[[1]]
inliers=hh[[2]]
print(paste0("Number of inliers: ",length(inliers)))
## [1] "Number of inliers: 11"
print("h=")
## [1] "h="
print(h)
##               [,1]          [,2]         [,3]
## [1,] -4.002669e-03 -8.798856e-04 -0.012050635
## [2,]  8.576245e-04 -3.954666e-03 -0.999896286
## [3,] -1.377376e-08  2.259075e-08 -0.005387763
par(mfrow=c(1,2))

plot(ima)
kpa %$% points(mx,my,col="red")
points(p1[inliers,],col="green")
plot(imb)
kpb %$% points(mx,my,col="red")
points(p2[inliers,],col="green")

Finally we apply the transformation to the first image and compare with the second one:

hm1=solve(h)

 map.affine <- function(x,y) {
    p=apply_homograph(hm1,cbind(x,y))
    list(x=p[,1],y=p[,2])
  }
  
nim=paste0(path,"set",set,"_",ida,".jpeg")
ima=imresize(as.cimg(readJPEG(nim)[,1:3099,]),scale=scl)
nim=paste0(path,"set",set,"_",idb,".jpeg")
imb=imresize(as.cimg(readJPEG(nim)[,1:3099,]),scale=scl)
  
imat=imwarp(ima,map=map.affine,dir="backward")

par(mfrow=c(1,2))
plot(imat)
plot(imb)

Now, after registration, we can compare the 2 images in order to detect differences (with a well defined time arrow if possible ;) )

d1=imat-imb
d2=(imat-imb)^2*(imat>0)
par(mfrow=c(1,2))
plot(d1)
plot(log(d2+.0001))

And that’s all for now. Probably the methodology can be fine tunned but is working prety well with most of the images. I hope this allow to the R coders to enter in the competition . Anyway, I think this workflow suggest some ways to attack the problem, build features or preprocess the images before feed them to a ML algorithm.

Good Image Registering!